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Abstract 

We study the calculation of the volume of the polytope B n of n x n 
doubly stochastic matrices; that is, the set of real non-negative matrices 
with all row and column sums equal to one. 

We describe two methods. The first involves a decomposition of the 
polytope into simplices. The second involves the enumeration of "magic 
squares" , i.e., nxn non- negative integer matrices whose rows and columns 
all sum to the same integer. 

We have used the first method to confirm the previously known values 
through n = 7. This method can also be used to compute the volumes of 
faces of B n . For example we have observed that the volume of a particular 
face of B n appears to be a product of Catalan numbers. 

We have used the second method to find the volume for n — 8, which 
we believe was not previously known. 

1 Introduction 

We study the calculation of the volume of the polytope B n of n x n doubly 
stochastic matrices; that is, the set of real nonnegative matrices with all row and 
column sums equal to one. This polytope is sometimes known as the Birkhoff 
polytope or the assignment polytope. We will describe and evaluate two meth- 
ods for computing the volume of B n . 

In the first method we decompose B n into a disjoint union of simplices all 
of the same volume and count the simplices. The fact that this can be done 



appears in [Stl]. This method applies to any face of B n as well. 

In the second method we count the number of n x n nonnegative integer ma- 
trices with all row and column sums equal to t (sometimes called magic squares) 
for suitable values of t. These numbers allow us to compute the Ehrhart polyno- 
mial of B n , which (essentially) has the volume of B n as its leading coefficient. It 
appears that this has been the most common method of computing the volume 
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of B n . Sturmfcls reports in iStu] on other work in which the volume of B n has 
been computed for n up to 7. We have also used this method to compute the 
volume when n = 8. 

This study is largely expository since the two methods are not new. However, 
the details about how we carry out these methods may be of interest. We are not 
aware of any reports of others who have carried out the simplicial decomposition 
method. 

As a byproduct of our program for carrying out the simplicial decomposition 
method, we are easily able to compute the volume of any face (of any dimension) 
of B n provided that n is not too large. This allowed us to discover that a certain 
special face of B n has a volume which appears to be given by a simple product 
formula. This formula is given in Conjecture |l|. 



Our study resulted from a question Mi] of Victor Miller, who asked how one 



could generate a doubly stochastic matrix uniformly at random. It is not hard 
to see that it would be easy to generate a random doubly stochastic matrix if 
one could easily calculate the volume of any face of B n . However the method 
described here for calculating face volumes is practical only for small n. 

In what follows we will make use of some well known properties of the face 
structure of B n : the vertices of B n are precisely the n\ n x n permutation 
matrices; on the other hand, for each pair (i,j) with 1 < i,j < n, the doubly 
stochastic matrices with entry equal to form a facet (maximal proper 





face) of B n and all facets arise in this way. See [BS] for further properties and 
references. 

In general, it is convenient to identify the faces of B n with certain n x n 
matrices of O's and l's, as follows. 

First we identify a 0-1 matrix with the set of entries in the matrix that are 
l's. Thus, for two 0-1 matrices A and B of the same size, we can define their 
union A U B as the 0-1 matrix whose set of l's is the union of the sets of l's of 
A and B. e.g., 



U 



Similarly we can speak of one 0-1 matrix containing another and so forth. 

Now to each face F of B n , we associate the matrix M which is the union 
of the vertices (permutation matrices) in F. The facets of B n containing F 
are precisely those associated with the zero entries of M. Since any face is the 
intersection of the facets containing it, any permutation matrix contained in M 
must be a vertex of F . Thus the vertices of F are precisely the permutation 
matrices contained in M, so we can recover F from M. In this way we identify 
the faces of B n with the set of 0-1 matrices which are unions of permutation 
matrices. Note that not every 0-1 matrix corresponds to a face of B n . For 
example 

1 

1 1 
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is not a union of permutation matrices, hence not a face of £?2- 



2 Volume 

It is easy to see that the dimension of B n is {n — l) 2 . Strictly speaking, the 
volume we wish to compute is the (n — l) 2 -volume of B n regarded as a subset of 
n 2 -dimensional Euclidean space. Thus, for example, the polytope Bi consists 
of the line segment joining the matrices 

1 \ / 1 

o i ) and { 1 

and hence its volume is 2. 

An n x n doubly stochastic matrix is determined by its upper left (n — 1) x 
(n — 1) submatrix. The set of (n — 1) x (n — 1) matrices obtained this way is 
the set A n of all nonnegative (n — 1) X (n — 1) matrices with row and column 
sums < 1 such that the sum of all the entries is at least n — 2. This is afHnely 
isomorphic to B n . In the Appendix we show that the ratio of the volume of B n 
to the volume of A n , regarded as a subset of Euclidean (n — l) 2 space, is n n ~ x . 
In some ways the volume of A n is easier to understand since its dimension is 
equal to the dimension of its ambient space. 



James Maiorana [Ma] (and probably others) noted a Monte Carlo method 
for approximating the volume of A n . Consider the set C n of (n — 1) X (n — 1) 
nonnegative matrices with row sums (but not necessarily column sums) < 1. 
This is the Cartesian product of n — 1 unit simplices in Euclidean (n — l)-space 
so its volume is ^ n _iy n -i ■ It is easy to choose points in C n uniformly at random. 
The probability a n that such a point is in A n is the ratio of the volume of A n 
to that of C n - Thus we can run Monte-Carlo trials to estimate a n and hence 
the volume of A n . 

For large n, this Monte Carlo method is impractical since a n is too small. 
However, it is useful for che cking computations for small n. A lower bound for 



a n is given by Bona in [Bo 



There is a more natural unit for the volume of B n and its faces. This is 
based on the fact that the vertices of B n are integer matrices. Suppose that F 
is a d-dimensional face of B n . Since its vertices have integer coordinates, the 
integer points in the affine span of F comprise a d-dimensional affinc lattice L. 
Given such a lattice there is a minimum volume of any d-simplex with vertices 
in L. Lattice points Wo, ... ,Wd are the vertices of one of these minimum volume 
simplices if and only if every point of L is uniquely expressible in the form 
Sf=o ki w ii where the k^s are integers whose sum is 1. The relative volume of a 
face F is the volume of F expressed in units equal to the volume of a minimal 
simplex in L. The relative volume of a face is the same whether regarded as a 
face of B n or as a face of A n , since the mapping from B n to A n (by taking the 
upper left (n— l)x(n— 1) minor) preserves integrality of points. 
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Here are the currently known relative volumes of B, 



n Relative Volume of B n 

1 1 

2 1 

3 3 

4 352 

5 4718075 

6 14666561365176 

7 17832560768358341943028 



8 12816077964079346687829905128694016 

To convert relative volumes to true volumes, we need to know the volume of 
a minimal simplex of A n . 

But the affine span of A n is all of (n — l) 2 -dimensional space. Hence the 
volume of a minimal simplex in A n is — i) a )i ' ana - ^he volume of a minimal 
simplex in B„ is 

3 Triangulations 

We call the first method for computing the volume of B n the triangulation 
method. The method applies to the calculation of the volume of any polytope. 
The essence is that we decompose the polytope into simplices and sum the 
volumes of the simplices. 

For B n we have used a standard method of decomposing a polytope P into 
simplices. See for example ptl| . To decompose P into simplices, we choose an 
arbitrary vertex v and form the collection of facets of P opposite v (facets of 
P not containing v.) We then recursively triangulate each facet. The triangu- 
lation of P is then formed by adding our chosen vertex to each simplex in the 
triangulation of each of the facets. 

The standard triangulations of B n and its faces have an unusual property, 
given in |Stl| , for which we provide a self-contained proof below. 

Proposition 1 In any standard triangulation of a face F of B n , every simplex 
has minimal volume in the affine lattice determined by F. 

Proof: Let F be a d-dimensional face of B n , vo any vertex in F, and G a 
facet of F opposite vq. Suppose that a simplex in a standard triangulation of G 
has vertices v±,V2, • ■ ■ , V4. We need to prove that the set of integer points of the 
affine space determined by F is the same as the set of points X^=o ^i v ii where 
the ki are integers whose sum is 1. 

Of course all the integer combinations are in the affine span. The question 
is whether there are any other points. 

Any integral point of the affine span can be uniquely expressed in the form 
Yli=o T i v ii where the r^'s are real numbers with sum 1. 
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Since vq is not in the face G, there is a facet of B n containing G but not Vo. 
Thus u must have at least one entry equal to 1 in the same position where all 
Vi, i > 1, have zeroes. Thus, in the hypothetical combination above, ro must be 
an integer. If we add ro(v\ — vo) to the combination above, we obtain another 
integral point in the affine span of G. It follows, using induction, that n + ro, 
and T2, ■ ■ ■ , rd are integers and therefore all the r's are integers, as desired. 

Note that, as a corollary, in any standard triangulation of a face of the B n , 
the number of simplices in the triangulation is equal to the relative volume. 

We also obtain an important computational principle. Given a face F of B n 
and a vertex v of F, the relative volume of F is the sum of the relative volumes 
of facets of F opposite v. 

4 The Triangulation Method for B n 

We now describe the triangulation method for computing the volume of B n . 
This is simply an elaboration of the principle that the relative volume of a face 
is the sum of the relative volumes of the faces opposite any vertex. 

We apply this principle recursively. To get started we use the fact that the 
relative volume of any zero-dimensional face of B n is 1. 

In the most naive plan we calculate the relative volumes of all faces. We first 
produce a list of all faces of each dimension. For dimension 0, we know all the 
relative volumes are 1. Then, for each face F of dimension d we select a vertex 
and find the opposite facets (of dimension d — 1). Assuming recursively that 
their relative volumes have already been computed, we now find the relative 
volume of F by summing the relative volumes of the facets. 

There are two serious drawbacks to the naive plan. 

Perhaps the most pressing problem is that we need to compute the volumes 
of an extremely large number of faces, since quite a few of the 2" possible 0-1 
matrices are actually faces of B n . Here we have recourse to a single important 
trick. If we permute the rows and columns of the matrix representing a face, we 
obtain the matrix of another face with the same volume. Also if we transpose 
a matrix representing a face, we obtain another face of the same volume. We 
regard matrices which can be obtained from each other by these operations as 
equivalent. We can cut down on the cost of our algorithm if we compute the 
volume only for a single "canonical" face in each equivalence class. 

The next most difficult problem is to produce the lists of faces. The most 
practical method that we found for producing faces is to start with the single 
(n — l) 2 -dimensional face, B n itself, and successively produce faces of lower 
dimension by intersecting with a facet of B n . While producing the faces we 
save the subface information so that we can look up the volumes when we are 
done. Unfortunately we need to construct a very large partially ordered set of 
faces before we can calculate any volumes since the only volumes we know are 
those of the zero-dimensional faces. While the cost in memory is not so bad 
for n less than 8, when we reach n = 8, we seem to need about 200 gigabytes 
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of intermediate storage. If the memory were available, the computation of the 
volumes would be relatively easy. In fact we are able to carry out a substantial 
fraction of the work before running out of memory. 

There are two phases to our algorithm. In the first phase we construct a 
collection of faces together with information about which ones are facets of which 
others. In particular, we successively compute, for d = (n— l) 2 , (n— l) 2 — 1, . . . , 0, 
a collection Td of c?-dimensional faces of B n . We begin by setting Ti n -iy = 
{B n }, i.e., consisting of just the all l's matrix representing B n itself. 

Given Td we produce Td-i as follows. Start with Td-i = 0- For each 
face F G Td we select a vertex v G F. We then find the facets of F opposite 
v, canonicalize these faces, and add them to Td-i- Having done this for all 
F G Td, we sort Td-i and remove the duplicates. Then, for each face / G Td-i, 
we save a list of pointers to the faces in Td from which / arose. (Equivalent 
faces can appear several times as opposite faces of the same face. When this 
happens, we include the pointer in the list of pointers multiple times.) 

This completes the first phase. In the second phase we start with To and 
work up to higher dimensions, calculating the relative volume of every saved 
face until we obtain the volume of B n itself. This is quite fast, requiring just 
one addition for each saved pointer. 

Note that once the pointers are constructed we do not need the faces them- 
selves, unless we want to know which face has each of the intermediate volumes 
we are computing. 

For larger values of n, the accumulators used for calculating the volumes 
will overflow. But we can get around this problem by using multiple precision 
arithmetic or by performing the volume calculation several times modulo various 
primes and combining the results with the Chinese Remainder Theorem. 

The main computational work of our algorithm takes place in three steps. 

1. for each face F G Td, find a vertex v G F. 

2. determine the facets of F opposite v. 

3. put these opposite facets into canonical form. 
We now describe how each of these steps is done. 

One important decision is the data structure for storing faces. We identify 
each face with the 0-1 matrix which is the union of its vertices (regarded as 
permutation matrices). For n < 8 it is convenient to represent each face as n 2 
bits of a single word, where the words of a (64-bit) computer are regarded as 
64- long arrays of bits. 

In Step 1 we are given a face / represented by a 0-1 matrix and we are 
looking for a permutation matrix 7r contained in /. This could be done with the 
assignment algorithm or one of the methods for finding maximum matchings, 
but for the small values of n that we were using, it was quicker to use a back- 
tracking search method, as follows. The matrix / has at least one 1 in its first 
row. We guess one of these as the location of the 1 in the first row of n. We 
then guess the location of the 1 in the second row of n, bearing in mind that it 
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cannot be in the same column as the 1 in the first row. We continue this way 
searching for the location of the 1 in subsequent rows. We backtrack if we reach 
a row in which there are no feasible choices. 

Now we consider Step 2. Given a face / and a vertex tt we need to find the 
facets of / opposite tt. 

For a moment let us ignore tt and consider the general problem of construct- 
ing the facets of /. The main principle is that each facet of / can be obtained 
by intersecting / with a facet of B n that does not contain /. Consider the facet 
corresponding to the pair The facet does not contain / if /y = 1. To in- 

tersect / with this facet we start by replacing with 0, obtaining a 0-1 matrix 
g. The face which is the intersection of / with the facet (i, j) is then the largest 
face h contained in g. The matrix h, which is the union of the permutation 
matrices contained in g, can be strictly contained in g. Given one of the l's in 
g, to test whether it is in h, we search for a permutation matrix in g which uses 
the 1 in question. This can be done with our backtracking search algorithm. 
The 1 in question is in h precisely when this search succeeds. When the position 
of a 1 in g is zero in h we say it is forced to zero. 

For example if n — 3 and / is the 3-dimensional face with matrix 



then the intersection of / with the facet corresponding to the middle entry of 
the top row is one-dimensional, with matrix 



In this example two zeroes in the last column are forced. This example also 
shows that although every facet of / is the intersection of / with a facet of B n , 
the converse is not true and the dimension of the intersection can be too small 
to be a facet of /. 

There are some additional simplifications when we search for the facets of 
/ opposite a given vertex tt of /. If g is a facet of / not containing tt, then g 
must contain a in place of one of the l's of tt. Thus there are at most n facets 
of / opposite tt. Observe that if g is a facet of / not containing tt, then g U tt 
is a union of permutation matrices and therefore a face of B n containing g and 
tt. Thus g U tt = f. This implies an important and helpful principle. When we 
introduce a at a 1 of tt and this results in a facet of /, then the only other 
positions that might be forced to zero are those of the other l's of tt. Thus we 
can loop through the n l's of tt one at a time and, for each of these, introduce a 
and determine what other 1 's of tt arc forced to be and produce accordingly 
a matrix, which we call a candidate. We obtain a set of n candidates among 
which all the facets opposite tt must occur. (This list can have duplicates which 
we remove.) Of these candidates the facets are those which are maximal under 
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inclusion. Indeed, it is clear that every candidate contains a face that has the 
same intersection with ir. But this face is contained in a facet which has an 
intersection with n that is at least as large. Thus every candidate is contained 
in at least one facet, and the facets are precisely the maximal candidates. 

Finally we describe Step 3, which we call "canonicalization" . 

The most straightforward way to choose a canonical form for a face / is to 
apply every element of our group of symmetries to / and choose the image of 
/ with the least value (where the bit pattern / is regarded as an integer.) But 
this is prohibitively slow. 

Instead we make use of certain special functions, which we call scores, which 
assign integers to every row and column of a 0-1 matrix. The scores have the 
special property that when rows are permuted, the row scores are permuted 
the same way leaving the column scores unchanged, whereas, when columns are 
permuted, the column scores are permuted the same way leaving the row scores 
unchanged. An example of an allowable score is to assign to each row its row 
sum and to each column its column sum. 

Given such scores we say that a matrix is in standard form if it satisfies the 
following three properties: 

1. the column scores are weakly increasing. 

2. the row scores are weakly increasing. 

3. in the case of tied row scores the rows are ordered lexicographically as bit 
strings. 

For a given 0-1 matrix, once its row and column scores have been computed 
it is easy to put a matrix and its transpose into standard form by forcing each 
of the three conditions above in the listed order. 

For each face constructed, we put both the face and its transpose into stan- 
dard form and finally choose the smaller of these two, regarded as integers, as 
the "canonical" form that is saved. 

Note that we are abusing terminology a little here since although the method 
always replaces a face by an equivalent face, it is conceivable that equivalent 
faces will canonicalize to distinct faces. When this happens, we still obtain 
correct volumes, but we end up doing work which could be avoided if the equiv- 
alence were recognized. However, if this event is rare, we obtain almost all the 
savings of true canonicalization as described above (but without the excessive 
cost). 

It turns out that just using row and column sums as the score functions fails 
to recognize a substantial number of equivalences. What we need are scores 
that tend to assign different values to different rows and columns. Slightly more 
complicated scores do better. Given a column score, we can produce a more 
complicated row score by assigning to each row the sum (or any symmetric 
function) of the values of the column scores of those columns for which l's 
occur in the given row. Similarly a row score can be used to produce a more 
complicated column score. We can also add two row scores to obtain another 



8 



row score or two column scores to obtain another column score. By combining 
steps like this we produced scores that were better at distinguishing rows (and 
columns) without being much more expensive to compute. 

This concludes our description of the triangulation method. As mentioned 
earlier it is reasonably practical for n < 8. The times required on a 500mhz 
DEC alpha were as follows: 



Time in seconds 
n < 6 less than 0.1 
n = 6 0.63 
n = 7 250.1 

Although the volumes of B n do not seem to follow a recognizable pattern, 
it seemed conceivable that there would be faces of B n for which the relative 
volumes had interesting properties. One fairly natural class is the set of matrices 
for which the set of of zeroes of the matrix form a Young tableau in a corner of 
the matrix. 

Since our triangulation method applies to any face of B n , we were able to 
check some natural classes of faces. It turned out that for the simplest non- 
trivial Young tableau faces the volumes apparently obey a simple rule, although 
we have not been able to supply a proof. More precisely, suppose that n > 2 
and that F n is the n x n matrix whose entry is 1 when j < i + 1 and 
otherwise. Then F n is a union of permutation matrices corresponding to a face 
of B n of dimension (™) with 2™ _1 vertices and we have the following 

Conjecture 1 The relative volume of F n is the product 




of the first n — 1 Catalan numbers. 



We have verified this for n < 12. 

Finally, we give some miscellaneous observations which may be useful but 
do not actually enter our algorithm. 

• In our method, we never needed to calculate the dimension of a face since 
the way they were produced guaranteed their dimension. However one 
may wonder how one can efficiently calculate the dimension of a face. 



One of the most efficient methods makes use of the fact, discussed in |BS|, 
that the dimension is equal to e+k — 2n, where e is the number of l's in the 
matrix of F and k the number of components in the graph corresponding 
to F. (i.e., the bipartite graph on 2n letters in which i is joined to j when 
the (i,j) entry of the matrix of F is 1.) 
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• The relative volume of any d-face F can be computed in several different 
ways since it is the sum of the relative volumes of the facets opposite any 
vertex of F. This yields linear relations on the volumes of (d — l)-faces of 
B n . It seems conceivable that these linear relations could be strong enough 
to yield useful information about the volumes. However from our limited 
investigation this does not appear to save anything in our computations. 

• Since our standard triangulations all involve minimum volume simplices, 
one might wonder whether all minimum volume simplices with vertices 
from the vertex set of B n belong to one of these triangulations. For n = 4, 
we found that there are 658584 minimum volume simplices whose vertices 
are vertices of B4. Of these, only 641112 belong to some standard trian- 
gulation. 



5 The Magic Squares Method 

In the next two sections we describe the magic squares method for calculating 
the volume of B n . We have no reason to believe that our implementation is 



subs tantially different from those used by others. (See |DG| ], |Mo| , | SS| , and 
[ (5tv| ] . ) The only apparent novelty is that we have carried out the computation 
when n = 8. 

We briefly explain here the connection between magic squares and the vol- 
ume of B n . 

It is known that for a rf-dimensional polytope P with integer vertices, for any 
nonnegative integer t, the number e(P, t) of lattice points contained in t ■ P is a 
polynomial of degree d in t. This polynomial is called the Ehrhart polynomial of 
P. Its leading coefficient is the volume of P in units equal to the volume of the 
fundamental domain of the affine lattice spanned by P. Thus if we know the 
values of e(P, t) for values of t from to d, we can find the Ehrhart polynomial 
by interpolation and in that way determine the volume of P. 

For P — B n , this method is particularly attractive since the polynomial is 
known to have certain symmetries, which make it necessary to calculate the 
values of e(B n ,t) for t only up to and including ( n ~ ) rather than (n — l) 2 . 

Note that e(B n , t) is exactly the number ofnxn matrices with nonnegative 
integer entries and all row and column sums equal to t, i.e., the number of n x n 
magic squares with sum t. In the next section we will describe how to count 
magic squares relatively efficiently. 

To see that we need only find e(B n ,t) for values of t up to and including 
("2 1 ) we re f er to the following identities: 

1. e(B n ,t) = for -n+1 < t < -1. 

2. e(B n , -n-t) = (-l)"- 1 ?^, t) for all t. 



These identities (conjectured in [ ADG|) are easy consequences of Ehrhart's Law 



of Reciprocity, which states that, for a c?-dimensional polytope P with integer 
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vertices, and t > 0, 

e*{P,t) = {-l) d e{P,-t) 

where e*(P, t) denotes the number of integer points in the interior of P. See |h| 
Chapter 9], and [|eJ for proof and references. 
Proof of 1 : 

e*(B n ,t) is the number ofnxn matrices with positive integer entries and all 
row and column sums equal to t. Since all the entries are > 1, each row and 
column sum must be > n, so e*(B n , t) = for 1 < t < n — 1. By Ehrhart's Law 
of Reciprocity this implies e(B n ,t) = for —n +1 < t < —1. 
Proof of 2: 

There is a one-to-one correspondence between n x n matrices with nonnegative 
integer entries and row and column sums t and n x n matrices with positive 
integer entries and row and column sums n + t. (Simply add 1 to each entry in 
matrices of the first type.) Thus e(B n ,t) = e*(B ni n + t). Applying Ehrhart's 
Law of Reciprocity, the right-hand-side equals (— e(B n , —n — t), which 
simplifies to (— l) n ~ 1 e(B n , — n — t). 

The effect of the first identity is that we know n — 1 zeroes of e(B n , t). We 
also have e(B n ,0) — 1. For each t > 0, if we calculate the value of e(B n ,t), 
by the second identity we obtain also the value of e(B n , —n — t). Thus if we 
calculate e(B n , t) for t up to (™), we have a total of n— l + l + 2(") = (n— 1) 2 + 1 
values of the e(B n ,t) so we have enough data to find the polynomial e(B n ,t) 
by interpolation. 



6 Counting Magic Squares 

We now describe the method we used for counting the number ofnxn magic 
squares of row and column sum t for t < ( n - 1 ). This seems no different from 



the methods used by others |DG| to carry out the smaller cases. 

Given an m-tuple r = (ri , . . . , r m ) and an n-tuple c = (ci , . . . , c n ) of nonneg- 
ative integers, we denote by N(r, c) the number of nonnegative integer matrices 
with row sums ri,...,r m and column sums ci, . . . , c„. 

There are a few computational principles. The first is that N (r, c) = 
unless |r| = J2i r i = J2j c j = l c l- Next note that N(r,c) is invariant under 
permutation of either the r's or the c's. Finally the principle that leads to 
substantial computational savings is that, for any integer k (usually near m/2) 

N(r, c)=J2 N((r u . . . , r k ), x)N((r k+1 ,. . . , 

X 

where the sum is over all nonnegative n-tuples x such that \x\ = r\ + • ■ • + r k 
and Xi < Ci, i = 1, . . . , n. This formula results from classifying the matrices 
counted by N(r, c) according to the column sums of the submatrix formed from 
the first k rows. For fixed column sums x\, . . . , x n , the column sums of the 
submatrix formed by the remaining rows must be Ci — Xi. The total number of 
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matrices in the class corresponding to x is the number of ways of choosing the 
top submatrix multiplied by the number of ways of choosing the bottom. 

The counting of magic squares amounts to the calculation of N(r, c) with 
the "constant" n-tuples r = c = (t, . . . , t). For this special case there are a few 
simplifications. We discuss the case when n is even. The same ideas apply with 
slight modification when n is odd. 

Suppose that n = 2m and we wish to calculate e(B n ,t). From our general 
principle we have 

e(B n ,t) = J2N(R,y)N(R,T-y) 

v 

where R is the m-tuple of all t's, T is the n-tuple of all i's, and y runs over 
all nonnegative n-tuples satisfying \y\ = mt, and yi < t for all i. For a fc-tuple 
y = (yi, . . . ,yk), let us denote by M(y) the number of distinct fc-tuples which 
arise by permuting the yi's. So, if z\,...,zi are distinct, and y is a fc-tuple 
consisting of k\ z\s, ki zi 's, etc., then M(y) — k\/(k\\ . . . fc;!). In terms of this 
notation a more computationally efficient version of the preceding equation is 

e(B n ,t)=J2M(y)N(R,y)N(R,T-y) (1) 
y 

where now we further restrict y to weakly increasing n-tuples. 

We can apply this principle again to the calculation of N(R, y) and N(R, T— 
y) that appear in the last formula. We find that 

N(R, y) = J2 M ( X ) N ( X > (Vu- ■■ , Vm))N(R - x, (y m+1 , . . . , y n )) (2) 

X 

where now x runs over all weakly increasing nonnegative m-tuples with |a;| = 
yi H h y m and x { < t for all i. 

We can save an additional factor of 2 by noting that the quantities N(R, y) 
are the same as N(R,T — y) except in a different order. Thus if we save the 
former in a suitable array, we can look up the latter ones in the array rather 
than computing them. 

Notice that the ingredients for calculating the sums N(R, y) and N(R, T — y) 
are the quantities N(x, y) where x and y vary over weakly increasing nonnegative 
m-tuples with Xi,yi < ("2 1 )- Thus it is sensible to precompute these quantities 
and save the results before forming the sums for N(R, y) or the sum for e(B n , t). 

For example, for n = 8 we need to precompute the quantities N(x, y) where 
x and y have length 4. Again it is easier to calculate 

N(x, y) = J2 N((x 1 ,x 2 ),z)N((x 2 ,x 3 ),y - z) 

where the sum is over all 4- long vectors z with \z\ = x\ + x 2 and Zi < yi 
for all i. However we do not have available the additional simplification to a 
sum over increasing sequences z. Thus on the right side we require the values 
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N(x, y), for pairs (x, y) where x has length 2 and y has length 4, not necessarily 
weakly increasing, where the components of x and y vary up to 21. It would 
be possible to precompute all the needed values and save these as well for later 
use. This might be advantageous since these results are used several times each. 
However, for simplicity, we use a subroutine to compute these, in effect repeating 
the calculation of any N(x, y) whenever needed. This subroutine in turn calls 
a subroutine for counting 2x2 matrices with prescribed row and column sums 
which calculates N ((x\, x 2 ), (2/1, 2/2)) = rnin (x 1, X2, yi, 2/2) + 1 whenever \x\ = \y\. 

The precalculation for n = 8 requires about 20 minutes on a 500Mhz DEC 
alpha. The remaining calculation also takes about 20 minutes. The first part 
can be calculated in single precision. In the remaining parts we need some 
sort of multiple precision method. We perform the calculation modulo several 
primes and combine the results with the Chinese Remainder Theorem. A similar 
program for n = 7 requires 38 seconds. 

Here are the Ehrhart polynomials e(B n , t), for n = 1, . . . , 8. For each n, the 
coefficient of the last binomial coefficient in the expression for e(B n ,t) is the 
relative volume of B n . We express the Ehrhart polynomial of B n as an integer 
combination of binomial coefficients C(t + n — 1 + k,n — 1 + 2k) = {^-l+tk)' 
k = 0, . . . , ("j 1 ) , because they are a basis for the polynomials satisfying p(— n — 

t) = (-lr-y*). 



e(Si,t) - C(t,0) 

e(B 2 ,t) = C(t + 1,1) 

e(B 3 ,t) = C(t + 2,2) + 3C(t + 3,4) 

e(B 4 ,t) = C(t + 3,3) + 20C(t + 4,5) + 152C(t + 5,7) + 352C(t + 6,9) 
e(B 5 ,t) = C(t + 4,4) + 115C(t + 5,6) + 5390C0 + 6,8) + 
101275C(t + 7, 10) + 858650C(t + 8, 12) + 
3309025C0 + 9, 14) + 4718075C(t + 10, 16) 
e(B 6 ,t) = C(t + 5,5) + 714C(t + 6,7) + 196677C(t + 7,9) + 
18941310C(t + 8, 11) + 809451144C(t + 9, 13) + 
17914693608C(t + 10, 15) + 223688514048C(t + 11, 17) + 
1633645276848C0 + 12, 19) + 6907466271384C(t + 13, 21) + 
15642484909560C(t + 14, 23) + 14666561365176C(t + 15, 25) 



13 



e{B 7 ,t) = C(t + 6,6) + 5033C(t + 7,8) + 

9090305C(i + 8, 10) + 4562637436C(t + 9, 12) + 
876755512997C(t+ 10, 14) + 80592643025748C(i+ 11, 16) + 
4085047594855542C(i + 12, 18) + 
125166504299043921C(t + 13,20) + 
2460507569635629206C(t+ 14, 22) + 
32199612314177385616C(t + 15, 24) + 
285953447105799237366C(t + 16, 26) + 
1 727929241 168643056768C(t + 17, 28) + 
6989369809320320632154C(t + 18, 30) + 
18096158896344747268932C(t+ 19,32) + 
27093648035077238674360C(t + 20, 34) + 
17832560768358341943028C(t + 21, 36) 

e(B s ,t) = C(t + 7,7) + 

40312C(i + 8,9) + 
544604804C(i + 9,ll) + 
1572522771472C(t + 10, 13) + 
1433860489078360C(i + 11, 15) + 
546197610013169408C(t + 12, 17) + 
104573799019751624800C(i + 13, 19) + 
11404657872578818785152C(t+ 14,21) + 
773100275338739807806336C(i + 15, 23) + 
34668602440014649185072000C(t + 16, 25) + 
1075823106306592550013512704C(t+ 17,27) + 
23865735845675030268755397632C(i + 18, 29) + 
387264682746696963082402212768C(t+ 19,31) + 
4666750907574155613393947915904C(t + 20, 33) + 
42107239094874587731729608526080C0+ 21, 35) + 
284859465667770778104594682157824C(i + 22, 37) + 
1435919936068954265096148477657088C(t + 23, 39) + 
5307981556350553774098942855517184C(t+ 24,41) + 
13958946247270195588626193027208192C(t+ 25,43) + 
24706461764218063045041689495950080C(t + 26, 45) + 
26368507913706408235698183181290240C(t + 27, 47) + 
12816077964079346687829905128694016C(t + 28, 49) 
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7 Comparison 



We now compare the two methods described above. 

The main advantage of the first method seems to be that it applies just as 
well to any face of B n as it does to B n itself. To apply the algorithm to a face 
F of B n , we simply start at the top level with the 0-1 matrix associated to F 
and produce lists of canonical subfaces as before. 

In the second method it is not obvious how well one could do in computing 
the volume of an arbitrary face F of B n . This would amount to counting the 
number of magic squares with prescribed zeros and row and column sums t for 
possibly as many as dim(F) values of t. We would not have e*(F, t) — for 
1 < t < n — 1, nor would we have e(F, t) = e*(F, n + t), because of the prescribed 
zeros in F. For certain F (e.g., those with the same number of prescribed zeros 
in every row) we would have an analogous identity, and some automatic roots, 
but in general we cannot guarantee any cutdown in the number of values of 
e(F, t) needed to determine the polynomial. Furthermore in the actual counting 
of magic squares with certain prescribed zeros, we would not be able to exploit 
the symmetries used in our algorithm above. 

The second method however has the advantage that, for computing volumes 
of B n itself, it is much more feasible in terms of memory. 

The second method also computes the Ehrhart polynomial. It seems possible 
that the first method could be modified to compute Ehrhart polynomials of the 
faces as well as just their volumes. We would need to keep track of the numbers 
of simplices of each dimension in a standard triangulation instead of just the 
simplices of the largest dimension. 
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Appendix: Ratio of Volumes of B n and A 



Consider the linear mapping C from (n — 1) x (n — 1) matrices to n x n matrices 
which sends matrix which is all zero except for a 1 at to the matrix 
which is all zero except for l's at and (n,n) and — l's at (i,n) and 

(n,j). If we follow £ by the addition of the n x n matrix that has the block 
form 

( J„_i,i \ 
y Ji,n-i 2 — n / 

where is the all fc x Z matrix of all l's, then we obtain the afBne mapping 
which sends A n to B n . Thus if we denote the ratio we seek by R, we find 
that R 2 is the determinant of the (n— l) 2 x (n — l) 2 matrix of dot products of 
F iiJ --F fciI =2 5 «*2*f. 

But in general if a;^ and are two m x m matrices and z is the m 2 x m 2 
tensor product matrix indexed by pairs ij and kl given by z^; = XikVji then 
detz = (detxdetj/) m . 

Our case is the special case that x = y = J„_i in _i + I n -i- Since the 
characteristic polynomial of — J m is A m_1 (A+m), the determinant of J„_i+7 n _i 
is n. It follows that R = ri n ~ 1 . 
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